NASA/CR- 1999-209361 
ICASE Report No. 99-27 



A High Order Discontinuous Galerkin Method for 2D 
Incompressible Flows 


Jian-Guo Liu 

University of Maryland, College Park, Maryland 
Chi-Wang Shu 

Brown University, Providence, Rhode Island 


Institute for Computer Applications in Science and Engineering 
NASA Langley Research Center 
Hampton, VA 


Operated by Universities Space Research Association 



National Aeronautics and 
Space Administration 

Prepared for Langley Research Center 
under Contract NAS 1-97046 


Langley Research Center 
Hampton, Virginia 23681-2199 


July 1999 



Available from the following: 


NASA Center for AeroSpace Information (CASI) 
7121 Standard Drive 
Hanover, MD 21076-1320 
(301)621-0390 


National Technical Information Service (NTIS) 
5285 Port Royal Road 
Springfield, VA 22161-2171 
(703) 487-4650 


A HIGH ORDER DISCONTINUOUS GALERKIN METHOD FOR 2D 
INCOMPRESSIBLE FLOWS 

JIAN-GUO LIU" AND CHI-WANG SHU+ 


Abstract. In this paper we introduce a high order discontinuous Galerkin method for two dimensional 
incompressible flow in vorticity streamfunction formulation. The momentum equation is treated explicitly, 
utilizing the efficiency of the discontinuous Galerkin method. The streamfunction is obtained by a standard 
Poisson solver using continuous finite* elements. There is a natural matching between these two finite element 
spaces, since the normal component of the velocity field is continuous across element boundaries. This allows 
for a correct upwinding gluing in the discontinuous Galerkin framework, while still maintaining total energy 
conservation with no numerical dissipation and total enstrophy stability. The method is suitable for inviscid 
or high Reynolds number flows. Optimal error estimates are proven and verified by numerical experiments. 
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1. Introduction and the Setup of the Scheme. We are interested in solving the following 2D time 
dependent incompressible Euler equations in vorticity streamfunction formulation: 

u>t + V • (uu?) — 0 

(1.1) Aft = u = 

u n — given on dQ, 

where V 1 - = ( — d y ,d x ). Notice that the boundary condition, plus the fact that u ♦ n = ^r, recovers ijy on 
the boundary (up to a constant) in a simple connected domain 

(1.2) V' |an= »/V 

We are also interested in solving the Navier-Stokes equations with high Reynolds numbers Re ^ 1: 

ujf H- V ■ (uad — ~ — Au 
Re 

(1.3) A = u;, u = 

u = given on <9fh 


The boundary condition is now (1.2) plus the non-slip type boundary condition: 


(1.4) 


dxi) 

dn 


= u 6 , r 


an 


For simplicity, we only consider the no-flow, no-slip boundary conditions v'b = 0, u^ r — 0 and periodic 
boundary conditions. 
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We first emphasize that, for Euler equations (1.1) and high Reynolds number (Re 3> 1) Navier-Stokes 
equations (1.3), it is advantageous to treat both the convective terms and the viscous terms explicitly. The 
methods discussed in this paper are stable under standard CFL conditions. Since the momentum equation 
(the tii st equation in (1.1) and (1.3)) is treated explicitly in the discontinuous Galerkin framework, there is 
no global mass matrix to invert , unlike conventional finite element methods. This makes the method highly 
efficient for parallel implementation, see for example [2]. As any finite element method, our approach has the 
flexibility for complicated geometry and boundary conditions. The method is adapted from the Runge- Kutta 
discontinuous Galerkin methods discussed by Cockburn et al. in a series of papers [7], [8], [9], [10], [11], [12], 
[20] and [6]. 

The main difficulties in solving incompressible flows are the incompressibility condition and boundary 
conditions. The incompressibility condition is global and is thus solved by the standard Poisson solver for 
the streamfunction </' using continuous finite elements. One advantage of our approach is that there is no 
matching conditions needed for the two finite element spaces for the vorticity u and for the streamfunction 
t/’. The incompressibility condition, represented by the streamfunction il>, is exactly satisfied pointwise, and 
is naturally matched with the convective terms in the momentum equation. The normal velocity u n is 
automatically continuous along any element boundary, allowing for correct upwinding for the convective 
terms and still maintaining a total energy conservation and total enstrophy stability. 

There is an easy proof for stability, both in the total enstrophy and in the total energy, which does not 
depend on the regularity of the exact solutions. For smooth solutions error estimates can be obtained. 

Our method, as it stands, can only compute 2D flows. Similar approach for the primitive variable 
formulation, suitable for 3D calculations, is under investigation. 

We do not advocate our method for modest or low Reynolds number flows. In such regime viscosity 
terms should be treated implicitly for efficiency. This is a much more challenging task in terms of space 
matching characterized by the Babuska-Brezzi-Ladyzenskaja condition, projection type methods, and global 
vorticity boundary conditions, see for example [3], [17], [25], [27], [18], [19], [24], etc. 

For convection dominated flows, as we are interested in this paper, we mention the work of Bell et al. 
[1] for second order Godunov type upwinding methods, see also Levy and Tadmor [22] and E and Shu [16]. 
This is still an active field for research. 

W’e now describe the setup of the scheme. We start with a triangulation T h of the domain D, consisting 
of polygons of maximum size (diameter) h, and the following two approximation spaces 

(L5 ) V* = {v-v |*€ P k (K), VA' € 7ft } , W k h = V h k n C 0 (fl), 

where P k (K) is the set of all polynomials of degree at most k on the cell K. 

For the Euler equations (1.1), the numerical method is defined as follows: find uj h € 1 ) k and V-’ft € Wk 
such that 

(1-6) & h v)k - u* * Vv) K + ^2 {u h • n CT h v~) e =0. Mv £ v£ , 

e£dK 

( 17 ) -(Vt/ift ■ Vy?) = (uh<p), V<p G , 

with the velocity field obtained from the stream function by 

(1-8) U ft = V x iph. 
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Here (•) is the usual integration over either the whole domain ft or a subdomain denoted by a subscript. 
Same thing for || ■ || for the L 2 norm. 

Notice that the normal velocity u/ ( n is continuous across any element boundary c, but both the solution 
LJh and the test function v are discontinuous there. We take the values of the test function from within the 
element K , denoted by v~ . The solution at the edge is taken as a single valued flux tU* , which can be either 
a central or a upwind biased averages For example, the central flux is defined by 


(1.9) u — 2 (f^h + ) 

where is the value of on the edge e from outside K , the complete upwind flux is defined bv 


( 1 . 10 ) 


LJf, — 


if u/j n > 0, 


^ uu'h if u h • n < 0. 

and the Lax-Friedriehs upwind biased flux is defined by 


( 1 . 11 ) 


Ufc • n Uj h = i [u/, • n + U h ) - a (u>,+ - u) h )] 


where a is the maximum of |u^ • n| either locally (local Lax-Friedrichs) or globally (global Lax-Friedrichs). 

We remark that, for general boundary conditions (1.2), the space 1 \ ^ h in (1.5) should be modified to 
take the boundary value into consideration. Moreover, additional physical vorticity boundary condition for 
any inlet should be known. 

Navier-Stokes equations (1.3) can be handled in a similar way, with the additional viscous terms treated 
by the local discontinuous Galerkin technique in [12], and with a local vorticity boundary condition in [13]. 
The detail is left to Sect. 3. Sect. 2 is devoted to the discussion of stability and error estimates for the Euler 
equations. Accuracy check and numerical examples are given in Sect. 4. Concluding remarks are given in 
Sect. 5. 

2. Stability and Error Estimates for the Euler Equations. For stability analysis, we take the 
test function v = ujf t in (1.6), obtaining 


J t \\ K||* - |<V • (wjju h )) K + 


Y k • nujh u) h ) e 

e€dK 


= 0, 


where we have used the exact incompressibility condition satisfied by u/, for the second term. Performing 
an integration by parts for the second term, we obtain 


d 1 
dt 2 


IK 111' + 


£<u h .n ( uj h uJ h 

e£dh 


- |K") 2 ))e = 0. 


Now, using the fact that 


U) = U) - (cj ) 2 — U ) 2 — uj[lo], 


where 


U) — ^ (u>~^ "f ^ ) , \<jj\ — — U) , 


we obtain 


d_ 1 
dt 2 


\\Wh\\ 2 K + ( Uh ' 11 {^h^h ~ ^l))e + \ ( Uh ,n Nm ” W A )) e = 0, 


etdK 


eedK 
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Notice that the second term is of opposite sign for adjacent elements sharing a common edge e, hence it 
becomes zero after summing over all the elements K (using the no-flow boundary condition on the physical 
boundary). The third term is the numerical dissipation: when Cft, is taken as the central flux (1.9), the third 
term is exactly zero; for the upwind flux (1.10), the third term becomes a positive quantity 

(2-D \ £<|u*.n|[u, h ] 2 > e 

tedK 

which is the total enstrophv dissipation. The effect of this is to control the size of the jump across the element 
interface and essentially “gluing” the solution there. Other upwind biased fluxes such as the Lax-Friedrichs 
flux (1.11) would produce a similar positive term as the total enstrophv dissipation. For smooth flows these 
jumps are of the order 0(h k ) within the truncation error of the scheme. We thus obtain the the following 
enstrophv inequality 

( 2 - 2 ) |lMI 2 ~ °’ 

which becomes an equality if the central flux (1.9) is used. 

The stability for the velocity field is now straightforward: we take p — %i> h in (1.7) to obtain: 

<W>A • W> k ) = - {UhVh) < ll*MIIK|| < C||V>/.||||a>fi|| 
by the Poincare inequality, which implies 

(2.3) IKH = HW'fcH < C\\uj h \\. 

Indeed, we can obtain a total energy conservation through the following arguments. Taking v = in (1.6), 
we obtain 


{d t u) h ip h ) h - (w/,Uft ■ V^/i)k + ^2 ( u h ■ n UJ h *l> h)e = 0, 

e edK 

Now the second term is zero since u* • — 0. The third term vanishes after summing over all elements 

since tl’h is continuous. Finally, noticing that 

-{d^n) = |l||V^|| 2 = |i||u ft || 2 , 

we obtain the conservation of energy 

( 2 - 4 ) | l | u „||=0 

even for a upwind flux. Thus there is no numerical dissipation for the energy. 

We now turn to the error estimates. For these we would need to assume that the solution is regular. 
Conceptionally, since this is a finite element method, the exact solution of the PDE satisfies the scheme 
exactly. As usual, we define the two projection operators: P is the standard I? projection into the space 
1 /* : and If is the standard projection into W* h : 

(V(v> - Ity) • Vy?) = 0, Vy>e IV* M . 

Denote the error functions by 

£ = U -U>h, S = ij) - %jj h 
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and their projections by 


Sfi — P £ — PjJ — UJh, 6f r — II 6 = II c - V/ ? . 


We first obtain a control of 6f } in terms of s: 

(VS-Vv) = -(e<p), Vv e Ho 1 ',, 

from the scheme (1.7) and the fact that the exact solution also satisfies (1.7). Now, taking ip = rf/,, we obtain 


(V6 h -V6 h ) = (X75-V6 h ) = ~(eS h ), 


which gives 


\m, I < C||e||. 


This leads to a bound for the velocity field 


llu - ufcii = II v(v» - ifo)ll < ||V(V; - m>)\\ + ||V(nv» - v*)ll 

(2.5) < ||V(V>-nvOII + C , ||e||. 

Since both the numerical solution and the exact solution satisfy (1.G), 

(2.6) {d t ev) K - ((wu - oj h u h ) ■ Vv) K + ^ ((u • nj - u,, • nCo, l )v~) e = 0. 


e€0K 


e v h k 


Take v — Sh . The second term becomes 

(2.7) ((wu - WhUft) • Ve h ) h - = (w(u - u h ) • V£ h )h + (eu ft • Ve) K - • V(u> - Pu>)) h :■ 

Noticing that u — u* is exactly divergence free, we may perform integration by parts to the first term on the 
right side of (2.7) to obtain 

(w(u - lift) • Ve h ) h - = (sk (u - Uft) ■ Vw) K + ((u - u h ) • nue^),.. 

eed A 

The second term on the right side of (2.7) is a complete derivative, hence can be integrated to give a pure 
boundary term 


(fTUfc • V£) K = i ( U h ' n ( £ ) 2 )r 


e£dK 


Plugging all these into (2.6) with v = e h , and collecting boundary terms, we obtain 
{d t £h£h)K + <£/,(u - u h ) • Vw) A - + (eu/,' ■ V(w - Pu>)) K + ^ I c 

eedK 


= 0 


where the boundary terms 

7 f = -((u - u h ) ■ n uJ£l) e - |((u h • n(£ _ ) 2 ) e + ((u-nu - u h ■ ntJ^) ep e 
= (u ft -n(££- - |(e~) 2 ))e 

= {Uft -n (££~ - i(0 2 ))p - i u h -nf (u> - (Poj)~))e . 
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Using the stability analysis in (2.2), we are left with 


(2.8) 


d 1 1 , . 

dt 2^ 


^E 


|-(e/i(u — u/, ) • Vu)) K - (s u/, • V(u? - Puj)) k 
+ E <u/,-n?(a;-(Pn;)-)) < .|. 


Assuming for the moment 


m IML<c. 

we ran first estimate* the boundary term 

E E <u fc • n?( W - (Pa-)-)),. < £ E C||[Pa;]|| e ||e1| e 

K eedh K teOK 

<M* + rX E iimi 2 - 

A eg 8K 

Using the above inequality together with (2.5) and (2.8), we now obtain 

|lM 2 < c ( Ikhll 2 + IIV(V^ - n^)|| 2 + ||a- - Puf m + ± E E IIMe) ■ 

V 1 A e£dK ) 

Here we understand the norms as a summation of the same norm on each K . Using the standard interpolation 
theory [5], we obtain 

J t \\Zh\\ 2 <C\\£ h f + Ch 2k 

which yields 


Ikftll < c h k . 


Together with (2.5), we have 

(2-10) ||u-Ufc|| + ||u;-w k || <Ch k . 

Using an inverse inequality, we have 


ll u - u fc||oo < Ch k 1 


this justifies the a priori assumption (2.9). 

The estimate (2.9) is optimal in terms of the space which is important since the main cost for the 
scheme is in the Poisson solver in W£ h . The vorticity estimate in (2.9) is however suboptimal with respect 
to the space U h *. If we use instead for the streamfunction and the upwind flux (1. 10), then a more 

detailed analysis will produce an order 0(h k+ 2 ) for the error in u, see [21] and [12] for details. However, we 
do not recommend this choice in practice, as the increase of half order accuracy is obtained with the price 
of one degree higher polynomials in the most expensive part of the algorithm, namely the Poisson solver. 
In our numerical experiments in Sect. 5, we observe that close to ( k + l)-th order of accuracy is generally 
achieved when fc-th degree polynomials are used in both the discontinuous space for and the continuous 
space for v, both for uniform and for non-uniform meshes. 
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3. The Scheme for the Navier-Stokes Equations. For the Navier-Stokes equations (1.3). there are 
two additional ingredients needing our attention: 

1. The viscous terms cannot be directly implemented in the discontinuous space VJ*. Instead, the stress 
tensor is first obtained locally using the same discontinuous Galerkin framework. 

2. Vorticity boundary values are not known physically. We obtain vorticity boundary conditions locally 
from the streamfunction using the kinematic relation in (1.3). 

We use the same finite element spaces and W { \ h defined in (1.5) for the vorticity and streamfunction, 
respectively. Denote, by the subspace of with zero value at the boundary nodes. Let Wjj be the 
finite element spaces extended from Wfi h with general non-zero values at the boundary nodes. The numerical 
method now becomes: 

(3.1) {d , ui h v) K - (w h u h ■ Vv)k + ^2 ( u /i ■ n ^i v ~)c = -{°h • Vv)k + (<7* • n v~) r , Vv € V£ h 

eedk eedK 

Notice that the test function is now in l^, (see [15]), and the stress tensor cr^ € (l^)^ is obtained from 
the vorticity bv the same discontinuous Galerkin framework: 

(3.2) Re (a h v) A = -(u) h V • v) A + ^ (u^v - • n) f , Vv € (Vjf )\ 

e£dK 

The fluxes and can be chosen as central averages 

(3-3) of, = k (/Th +fT h)' w* = 5 ( w /7 + -7 ) 

or better still, as alternate one-sided fluxes, namely, at each edge e with an arbitrarily fixed orientation, one 
of Oh and uJ/* is taken as the left value and the other taken as the right value. It can be verified that, for 
k — 0 and a rectangular triangulation, the central fluxes (3.3) produce a wide stencil central approximation 
to the second derivatives (uJi- 2 , & i and uij +2 are used for approximating cj^), while the alternate one-sided 
fluxes produce a compact stencil central approximation (cj,_i, ujj and i are used for approximating o; w ). 
Also, numerical and theoretical evidence shows that the alternate one-sided fluxes produce more accurate 
results [12]. In this paper we use only the alternate one-sided fluxes for the viscous terms. 

We advocate the same steps as in [15] for a finite element method. In the time stepping, we first compute 
the vorticity at the interior nodes, and we will use these values to compute a stream function and then we use 
the stream function to determine the vorticity at the boundary nodes. This time-stepping is very efficient 
and we do not need any iteration between the boundary vorticity value and interior values, thus eliminating 
some traditional difficulties associated with the vorticity formulation. This time-stepping was first developed 
for finite-differences in [13, 14]. 

Since (3.1) is treated explicitly, The value of u;" +1 is computed via two steps. First, we compute 
(cj£ +1 v)k for all v 6 V^ h from the explicitly time stepping of (3.1). The value of u;£ +1 at the interior 
element can now be directly computed from this term. However, the value of at the boundary element 
shall be determined after we computed the stream function as we explain below. 

These values, (a;£ +1 v)k for all v € is sufficient to be used to compute the stream function from 

(3.4) - W +1 • V<f) = K +1 V ), V<p € W * h , 
with the velocity field obtained from the stream function by 

(3.5) U jf +I =V- L V£ +l - 
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We now describe how to get the vorticity at the boundary cells. Since is known, we can compute 

-<v< +1 • v v ) = K +1 *>), 

for the* test func tion ^ at the boundary nodes. can then use it to compute the value of vorticity at the 
boundary elements. 

For problems with periodic boundary conditions, the formulation above admits the following stability 
results: 

(3-6) ^M 2 + 2|K|| < 0, 

which in turn implies stability for the velocity field (2.3). The proof is similar to the Euler case, see [12] 
for the details. W ith the vorticity boundary condition mentioned above, we are unable to obtain a stability 
estimate. However, this type of vorticity boundary treatment for conventional finite difference and finite 
elements is stable, see [14] and [15]. 

4. Accuracy Check and Numerical Examples. W ; e implement our method on triangulations based 
on rectangles. W hen a P k result is referred to it is obtained with P k elements for the? vorticity u) and Q k 
elements for the streamfunction '</•’. Strictly speaking Q k elements should also be used for the vorticity u 
for the exact energy conservation (2.4) to hold, however to save cost we use P k elements for the vorticity 
^ instead. Energy stability (2.3) and enstrophy stability (2.2) still hold in this case. We have used both 
the upwind flux (1.10) and the (global) Lax-Friedrichs flux (1.11) for the calculations, however we will only 
show the results obtained with the Lax-Friedrichs flux to save space. The time discretization is by the third 
order positive Runge-Kutta methods in [26]. 

Example 1 : This example is used to check the accuracy of our schemes, both for the Euler equations (1.1) 
and for the Navier-Stokes equations (1.3) with Re — 100, for both the periodic and the Dirichlet boundary 
conditions, and with both a uniform mesh and a non-uniform mesh. The non-uniform mesh is obtained by 
alternating between 0.9 Ax and 1.1 Ax for the mesh sizes in the x direction, similarly for the mesh sizes in 
the y direction. The initial condition is taken as 

(4-1) w(x.y.O) - —2 sin(x) sin(y), 

which was used in [4]. The exact solution for this case is known: 

(4.2) tj(x, y,t) = — 2sin(x)sin (y)e~&. 

We use the domain [0, 2n] x [0,27r] for the periodic case and [0, tt] x [0,7r] for the Dirichlet case and compute 
the errors at t = 2 for the periodic case and at t = 1 for the Dirichlet case. We list in Table 4.1 (uniform 
mesh) and Table 4.2 (non-uniform mesh) the L\ and L ^ errors, at t = 2, measured at the center of the cells, 
for the periodic boundary conditions. Table 4.3 (uniform mesh) and Table 4.4 (non-uniform mesh) contain 
the results with the Dirichlet boundary conditions at t = 1. W 7 e remark that, because of the difference in 
the sizes of the domains of the periodic and Dirichlet cases, the errors with the same number of cells are of 
different values, but the orders of accuracy are similar. W 7 e have also computed the errors of the relevant 
derivatives at the centers of the cells, which help in giving us truly Loc errors throughout the domain. We 
will not show them to save space. 
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Table 4.1 

Accuracy test, uniform meshes, periodic boundary conditions. 



Euler 

Navier-Stokes with Re = 100 

mesh 

L { error 

order 

L 00 error 

order 

L l error 

order 

L 00 error 

order 

PI 

16 2 

7.77E-03 


1.80E-02 


7.65E-03 

— 

1.82E-02 


32 2 

1.01E-03 

2.94 

2.46E-03 

2.87 

1.03E-03 

2.89 

2.55E-03 

2.83 

64 2 

1.28E-04 

2.99 

3.14E-04 

2.97 

1.36E-04 

2.92 

3.44E-04 

2.89 

128 2 

1.60E-05 

3.00 

3.94E-05 

2.99 

1.80E-05 

2.92 

4.63E-05 

2.89 

P2 

16 2 

6.26E-04 

- 

1.58E-03 

-- 

2.0GE-04 

— 

5.85E-04 

— 

32 2 

5.52E-05 

3.50 

2.75E-04 

2.52 

1.37E-05 

3.90 

3.24E-05 

4.17 

64 2 

4.82E-06 

3.52 

3.8 IE-05 

2.85 

2.40E-0G 

2.51 

4.10E-06 

2.98 

128 2 

4.04E-07 

3.58 

4.96E-06 

2.94 

4.05E-07 

2.57 

6.44E-07 

2.67 

P3 

16 2 

9.74E-05 

— 

2. 3 IE-04 

— 

9.68E-05 

— 

2.33E-04 

- - 

32 2 

6.81E-06 

3.84 

1.67E-05 

3.79 

6.22E-06 

3.96 

1.50E-05 

3.9G 

64 2 

4.36E-07 

3.96 

1.05E-06 

3.99 

3.82E-07 

4.02 

9.25E-07 

4.02 

128 2 

2.71E-08 

4.01 

6.59E-08 

3.99 

2.33E-08 

4.04 

5.70E-08 

4.02 


Table 4.2 

Accuracy test , non-untform meshes , periodic boundary conditions. 



Euler 

Navier-Stokes with Re = 100 

mesh 

L l error 

order 

L°° error 

order 

L 1 error 

order 

L°° error 

order 

pi 

16 2 

8.49E-03 

— 

2.85E-02 

— 

7.77E-03 


2.80E-02 

_... 

32 2 

1.44E-03 

2.5G 

5.56E-03 

2.36 

1.16E-03 

2.75 

5.45E-03 

2.36 

64 2 

2.81E*04 

2.3G 

1.13E-03 

2.29 

2.17E-04 

2.42 

1.03E-03 

2.40 

128 2 

5.90E-05 

2.25 

2.59E-04 

2.13 

4.13E-05 

2.40 

1.94E-04 


P2 

mm 

7.88E-04 

— 


— 

3.37E-04 

— 

1.18E-03 

— 


7.82E-05 




1.78E-05 


6.40E-05 


64 2 

7.66E-06 

Kjfi 







CM 

00 

CM 

r-H 

7.43E-07 

3.37 

6.11E-06 

3.07 

4.34E-07 

2.60 

1.00E-06 

2.80 

P3 





— 


— 


— 




2.60E-05 

3.65 


3.96 


3.92 

64 2 

4.60E-07 

3.9G 

1.77E-06 

3.88 

4.01E-07 


1.28E-06 

4.06 

128 2 

2.86E-08 

4.01 

1.09E-07 

4.02 

2.44E-08 

4.03 

7.70E-08 

4.06 


9 














Table 4.3 

Accuracy test , uniform meshes , Dirichlet boundary conditions. 



Euler 

Navier-Stokes with Re = 100 

mesh 

L 1 error 

order 

L x error 

order 

L l error 

order 

L°° error 

order 

pi 

n 


IKS 


— 

5.75E-04 

— 

1.32E-03 

— 


8.19E-05 

2.85 

1.92E-04 

2.68 



1.78E-04 

2.89 

64- 

1.06E-05 



mm 


mam 

3.76E-05 

2.25 


1.35E-06 

2.98 

1.42E-05 

1.92 

1.25E-06 

2.94 

8.06 E-06 


P2 

16- 

4.76E-05 

— 

2.57E-04 

— 

1.51E-05 

-• 

4.05E-05 

— 

32 L> 

4.28E-06 



2.85 

2.49E-06 

2.60 

6.09E-06 

2.73 

64 2 




2.94 

4.11E-07 

2.60 

9.42E-07 

2.69 

128 2 

3.17E-08 

3.5G 

5.92E-07 

2.97 

6.16E-08 

2.74 

1.34E-07 

2.81 

P3 

16-’ 

6.80E-06 


1.58E-05 

— 

6.34E-06 

— 

1.53E-05 

— 

32 2 

4.22E-07 

4.01 

1.06E-06 

3.90 

3.90E-07 

4.02 

9.45E-07 

4.02 

64 2 

2.66E-08 

3.99 

6.90E-08 

3.94 

2.38E-08 

4.04 

5.81E-08 

4.02 

128 2 

1.G6E-09 

4.00 

4.25E-09 

4.02 

1.46E-09 

4.03 

3.59E-09 

4.02 


Table 4.4 

Accuracy test, non-unifonn meshes , Dirichlet boundary conditions. 



Euler 

Navier-Stokes with Re — 100 

mesh 

L l error 

order 

L°° error 

order 

L} error 

order 

L°° error 

order 

pi 

16 2 

1.12E-03 

— 

4.35E-03 

— 

9.93E-04 

— 

4.25E-03 


32 2 

2.44E-04 

2.20 

9.79E-04 

2.15 

1.95E-04 

2.35 

8.74E-04 

2.28 

64 2 

5.61E-05 

2.12 

2.39E-04 

2.04 

3.90E-05 

2.33 

1.75E-04 

2.32 

128 2 

1.36E-05 

2.04 

6.29E-05 

1.92 

7.76E-06 

2.33 

3.54E-05 

2.31 

P2 

16 2 

7.54E-05 

N9 


— 

1.98E-05 

-- 

6.63E-05 

.... 

32 2 

8.15E-06 

mu 

4.33E-05 

2.93 

2.61E-06 

2.93 

7.03E-06 

3.24 

64 2 

8.46E-07 

3.27 

5.35E-06 




1.02E-06 

■ESll 

128- 

8.31E-08 






1.49E-07 

2.78 

P3 

16 2 

7.17E-06 

— 

2.49E-05 

— 

6.65E-06 

— 

2.18E-05 

— 

32 2 

4.46E-07 

4.01 

1.66E-06 

3.91 

4.09E-07 

4.02 

1.31 E-06 


64 2 


3.99 

1.04E-07 

3.99 

2.50E-08 

4.03 

7.86E-08 

Ki!Sil 

CN 

00 

CN 

1.75E-09 

4.00 

6.89E-09 

3.92 

1.53E-09 

4.02 

4.77E-09 

4.04 
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We can clearly see from these tables that dose to (fc-f- l)-th order of accuracy is generally achieved when 
k - th degree polynomials are used in both the discontinuous space' for u) and for the Poisson solver, both for 
the uniform and for the non-uniform meshes. 


Example 2 : The double shear layer problem taken from [1]. We solve the Euler equation (1.1) in the 
domain [0, 2n] x [0. 2 tt] with a periodic boundary condition and an initial condition: 


Scos(t) - i serh 2 ((y - n/2)/p) y < tt 
<5 cos(x) + i se.ch 2 ((3ir /2 — y)/f>) y > tt 


where we take p = 7t/ 15 and 6 — 0.05. 

The solution quickly develops into roll-ups with smaller and smaller scales, so on any fixed grid the full 
resolution is lost eventually. We use fixed uniform meshes of 64 x 64 and 128 x 128 rectangles and perform 
the calculation up to t = 8. We plot the time history of total energy (square of the L 2 norm of velocity u) 
and total enstrophy (square of the L 2 norm of vorticity a;) in Fig. 4.1, as well as contours of the vorticity jj at 
t = 6 in Fig. 4.2 and at t = 8 in Fig. 4.3 to show the resolution. We can see from Fig. 4.1 that the numerical 
dissipation decreases roughly in the order of P l 64 2 , P l 128 2 , P 2 64 2 , P 3 64 2 , P 2 128 2 , and P 3 128 2 . The 
higher order methods have better resolutions and in general the resolution is quite good judging from the 
contours. We remark that when the numerical viscosity becomes too small with higher order methods, since 
the schemes are linear, numerical oscillations are unavoidable when resolution to sharp fronts is lost., leading 
to instability. This is common for all linear schemes. However, the discontinuous Galerkin method we use 
here is able to get stable solutions for much sharper fronts with the same mesh than central type finite 
difference or finite element methods. More extensive numerical resolution study for this example can be 
found in [23]. For a comparison with nonlinear ENO schemes, we refer to [16]. 




Fig. 4.1. The time history of energy (square of the L 2 norm of the velocity u) and total enstrophy (squai'e of the L 2 norm 
of vorticity u) ) . P 1 with 64 2 mesh in solid line, P 1 with 128 2 mesh in dashed line , P 2 with 64 2 mesh in dash-dot line, P 2 with 
128 2 mesh in dotted line, P 3 with 64 2 mesh in long dashed line , and P 3 with 128 2 mesh in dash-dot-dot line . 


Example 3 : The vortex patch problem. We solve the Euler equation (1.1) in [0,27 t] x [0,27t] with the 


ll 



following initial condition: 


(4.4) 


uj( x, y. 0 ) 


' -1- f <*< 'f, f <y< 

< h f < * < t 1 , 

0, otherwise 


and periodic boundary conditions. The contour plots of vorticity u,\ with 30 equally spaced contour lines 
between = -1.1 and u? = 1.1, are given in Fig. 4.4 for t, = 5 and in Fig. 4.5 for t = 10. We can see that 
the scheme gives stable results for all runs, and higher order schemes give better resolutions for vorticity. 

5. Concluding Remarks. We have developed a high order discontinuous Galerkin method for the two 
dimensional incompressible Euler and Navier-Stokes equations in the vorticity streamfunction formulation, 
coupled with a standard continuous finite element solution of the Poisson equation for the streamfunction. 
A natural matching between the two finite element spaces allows us to obtain total energy conservation and 
total enstrophy stability. Numerical examples are shown to demonstrate the accuracy and resolution of the 
methods. 


Acknowledgments. This work was initialized when both authors were visiting the Mittag-Leffier In- 
stitute in Sweden. We thank the faculty and staff there for their warm hospitality. 
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Fig. 4.2. Contour of vorticity u ) at t — 6. 30 equally spaced contour lines between = —4.9 and lj = 4.9. 
with 64 2 mesh; Right: results with 1 28 2 mesh. Top: P l ; middle: P 2 , bottom: P*. 


Left: results 
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PI, 64* mash, t=6 

30 equally spaced vofdclty contour* from -4.0 to 4.0 


PI, 120* mart, t=S 

30 aq unity spa cad vort lefty contour* from -4.0 to 4.0 



Fig. 4.3. Contour of vorticity uj at t — 8. 30 equally spaced contour lines between lo = —4.9 and uj = 4.9. Left: results 
with 64 2 mesh; Right: results with 1‘28 2 mesh. Top: P 1 ; middle: P 2 , bottom: P 3 . 
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Vorticity at t=5, PI , 64 2 mesh 

30 equally spaced contours from -1 .1 to 1 .1 
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Vorticity at t=5, P2, 64 2 mesh 

30 equally spaced contours from -1 .1 to 1 .1 
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! Vorticity at t=5, P3, 64 2 mesh 
j 30 equally spaced contours from -1 .1 to 1 .1 





Fig. 4.4. Contour of vorticity u) at t ~ 5. 30 equally 
with 64 2 mesh; Right: results with 128 2 mesh. Top: P 1 ; m 


Vorticity at t=5, PI , 1 28 2 mesh 

30 equally spaced contours from -1.1 to 1 .1 



iced contour lines between uj = - 1.1 and oj = 1 . 1 . left: results 
lie: P 2 , bottom: P 3 . 
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Vorticity at t=1 0, PI , 64 2 mesh 
30 equally spaced contours from -1 .1 to 1 .1 
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Vorticity at t=1 0, P3, 64 2 mesh 

30 equally spaced contours from -1 .1 to 1 .1 
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Vorticity at t=1 0, PI , 1 28 2 mesh 
30 equally spaced contours from -1 .1 to 1 .1 
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Vorticity at t=10, P2, 128 2 mesh 
30 equally spaced contours from -1 .1 to 1 .1 



Vorticity at t=1 0, P3, 1 28 2 mesh 
30 equally spaced contours from -1.1 to 1 .1 



1 1 i i 1 1 1 1 1 1 1 1 1 1 1 1 


Fig. 4.5. Contour of vorticity u : at t — 10. 30 equally spaced contour lines between u> = —1.1 and uj — 1.1. Left: results 
with 64 2 mesh; Right: results with 128 2 mesh. Top: P 1 ; middle: P 2 , bottom: P 3 . 
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